# 找到与表型相关的模块

# Rscript step5.R --exprFile "../test_data/datExpr.csv" --traitsfile "../test_data/datTraits.csv" --network "../1_construct_network/network.RData" --traitName weight_g --output "../5_find_key_drivers"

library(optparse)


option_list = list(
  make_option("--exprFile", type="character", default=NULL,
              help="expression file path"),
  make_option("--traitFile", type="character", default=NULL,
              help="traits file path"),
  make_option("--network", type="character", default=NULL,
              help="network file path"),
  make_option("--traitName", type="character", default=NULL,
              help="trait name"),
  make_option("--output", type="character", default="out.txt",
              help="output path [default= %default]")
)

# 解析命令行参数
opt_parser = OptionParser(option_list=option_list, add_help_option=TRUE)
opts = parse_args(opt_parser)

datExpr = read.csv(opts$exprFile, stringsAsFactors = FALSE, row.names = 1);
datTraits = read.csv(opts$traitFile, stringsAsFactors = FALSE, row.names = 1);
traitName = opts$traitName

load(opts$network)

library(WGCNA)

nGenes = length(net$goodGenes)
nSamples= length(net$goodSamples)
moduleColors = labels2colors(net$colors)

# Define variable weight containing the weight column of datTrait
weight = as.data.frame(datTraits[[traitName]]);
names(weight) = "weight"


MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs = orderMEs(MEs0)

modNames = substring(names(MEs), 3)
geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p"));
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples))

names(geneModuleMembership) = paste("MM", modNames, sep="");
names(MMPvalue) = paste("p.MM", modNames, sep="");
geneTraitSignificance = as.data.frame(cor(datExpr, weight, use = "p"));
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));
names(geneTraitSignificance) = paste("GS.", names(weight), sep="");
names(GSPvalue) = paste("p.GS.", names(weight), sep="");


GS = abs(geneTraitSignificance$GS.weight)
png(
  filename = paste(opts$output, paste(traitName, "_modules_significance.png"), sep="/"),
  width = 18,
  height = 9,
  units = "in",
  res = 300
)
par(mfrow=c(1,2))
cex1=0.9

plotModuleSignificance(GS, moduleColors)

# 求均值
moduleGS <- tapply(GS, moduleColors, mean)
module_significance <- data.frame(
  Module = names(moduleGS),
  Mean_GS = as.numeric(moduleGS),
  Size = as.numeric(tapply(GS, moduleColors, length))
)

# 按平均 GS 值排序
module_significance <- module_significance[order(module_significance$Mean_GS, decreasing = TRUE), ]
max_sign_module = module_significance$Module[1]

print(max_sign_module)


column = match(max_sign_module, modNames);
moduleGenes = moduleColors==max_sign_module;

verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
                   abs(geneTraitSignificance[moduleGenes, 1]),
                   xlab = paste("Module Membership(", max_sign_module, ")"),
                   ylab = paste("Gene significance(", traitName, ")"),
                   main = paste("Module membership vs. gene significance\n"),
                   cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = max_sign_module)


dev.off()

# save result
write.csv(module_significance,
          file = paste(opts$output, paste(traitName, "_modules_significance.csv"), sep="/"),
          row.names = FALSE)
